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Planar Solid-Oxide Fuel Cell Research and Development 

S. Kapadia, W. K. Anderson, J. Newman, B. Hilbert 

University of Tennessee SimCenter at Chattanooga, 

701, East M.L. King Boulevard, Chattanooga, TN - 37403, U.S.A. 

Abstract 

Numerical simulation capabilities to perform flowfield analysis and sensitivity 
analysis of solid oxide fuel cells (SOFCs) have been developed. A three-dimensional, 
implicit, multi-species fuel cell solver is modified to analyze cell geometries supplied by 
NexTech Materials (NexTech). Four different NexTech cell configurations of varying 
complexities have been analyzed in this study. A capability to perform thermo-mechanical 
analysis of NexTech cells is developed by coupling the fuel cell code with the structures 
code. The capability was initially developed for a simplified cell that also included 
sensitivity analysis of thermo-mechanical variables. Sensitivity derivatives are computed 
using the direct differentiation method in the fuel cell and structures codes. 

1. Introduction 

Development of alternative energy producing devices has attracted a lot of attention 
from researchers all over the world in recent years. Due to its capability of producing 
efficient energy using various hydrocarbon fuel types. Solid Oxide Fuel Cells (SOFCs) 
present a promising technology for the future. Even though SOFCs are still in the 
developmental stage, numerical techniques can be effectively utilized to find solutions to 
many of the design hurdles affecting the commercial application of the technology. 
Experimental and numerical approaches have been undertaken by several researchers [1-8] 
to study the behavior of SOFC. Some of the advantages of the numerical simulations over 
the experiments arc the cost effectiveness and the fact that the simulations provide a wealth 
of data that is difficult or impossible to obtain experimentally and can be used to perform 
in-depth analysis of the SOFC unit/system. Numerical simulations can contribute greatly 
toward the development of better designs that can produce more power, increased 
efficiency, and extended life expectancy of various SOFC components. 

To date, numerical simulations have been primarily focused on analysis of fuel 
cells or fuel cell components, without strong emphasis on utilizing the simulations in a 
design optimization environment. Because of the emphasis on analysis instead of design, 
sensitivity information to determine the effects of variations in design parameters on 
performance has been primarily implemented by simply changing the parameter of interest, 
re-running the simulation, and comparing the results with those from the original 
simulation [1,5-6,8]. While this approach can be used to determine the effects of 
parameter variations on fuel cell performance, a more rigorous approach toward 



optimization would likely lead to better designs, and can also provide improved insight 
into the parameters affecting the performance of the fuel cell. For SOFC problems, 
example cost functions that can be used for improving performance include minimizing 
temperature variations, obtaining equal distribution of fuel in each of the channels, 
minimizing stress inside different components or maximizing power. Design variables 
may be related to the shape/size of the fuel channels, electrodes, electrolyte, and 
interconnect, but may also be coupled to the stoichiometric composition of fuel or material 
properties such as the porosity or tortuosity of the electrodes. 

In references [9] and [10], optimization algorithms have been used to improve the 
performance of a polymer-electrolyte-membrane fuel cell (PEM) using four design 
variables, where the sensitivity derivatives used for the optimization algorithm have been 
obtained using a finite-difference approach. While finite differences are often a viable 
means for computing sensitivity derivatives, this method can be computationally restrictive 
when a sufficiently large number of design variables are present. In addition, accurate 
derivatives can sometimes be difficult to obtain using finite differences because of 
subtractive cancellation errors [11], which occur when the function evaluations in the 
numerator become computationally indistinguishable [12] when very small perturbations 
are used. By using a direct differentiation or discrete adjoint method [16-27], sensitivity 
derivatives that are consistent with the flow solver may be obtained for use in a design 
optimization environment. The fuel cell code utilized in this report is capable of 
computing sensitivity derivatives of a cost function with respect to desired parameters 
using both direct differentiation and discrete adjoint methods. Although sensitivity studies 
shown in this report uses direct differentiation method to compute sensitivity derivatives, 
implementation and verification of discrete adjoint method in the fuel cell code has been 
published in previous work [28-31]. In different design studies performed at the 
SimCenter, sensitivity derivatives were also utilized in a design environment to optimize 
cost functions representing uniform fuel distribution [28,31], temperature distribution [29] 
and cell voltage [31]. 

Due to the presence of various thermal mechanisms inside SOFC, temperature 
distribution exhibits non-uniformity throughout the domain. This coupled with the 
mismatch in coefficients of thermal expansion of different SOFC components makes stress 
analysis an interesting avenue to explore. Recently, few studies have been performed [32- 
38] to analyze stress components inside different components of SOFC. Lin et al. [32] 
analyzed effects of clamping load on the thermal stress distribution in a planar SOFC. Five 
different compressible loads were applied to investigate effects on stress distribution. 
Gulfam et al. [35] analyzed thermal stress inside PEN region of SOFC for co-flow, 
counter-flow and cross-flow configurations using commercial software, ABAQUS [39]. 
Maximum stress inside anode layer was found at high-temperature regions located on the 
anode-electrolyte interface for all configurations. Jiang et al. [36] performed thermal stress 
analysis of SOFC with the bonded compliant seal design. Stress analysis performed using 
commercial software, FLUENT [40] and ANSYS [41] investigated effects of temperature 
non-uniformity and cell voltage on thermal stress distribution. Weil and Koeppel [37] also 
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analyzed effects of different seal designs on stress-strain state. They used glass-ceramic, 
brazed joints and foil-based seals in this study. Chiang et al. [38] analyzed effects of anode 
porosity on thermal stress in anode-supported SOFC using commercial software, STAR- 
CD [42] and MARC [43]. The study indicated presence of higher principal stress at low 
cell voltages due to high local current density and steep temperature gradients. 

None of the studies mentioned in the previous paragraph attempted to compute 
sensitivity derivatives of a cost function involving thermal stress components using formal 
procedures such as direct differentiation or discrete adjoint methods. Such capability is 
required if optimization of principal stress components or strain rates with respect to 
geometrical or material parameters is desired. One of the cases (Case 5) described in this 
report demonstrates such capability, where stress sensitivities are computed with respect to 
the cathode porosity. Implementation of this method is not trivial, as it requires 
computation of sensitivity derivatives of flowfield variables in the fuel cell code and 
coupling them with the structures solver to compute stress sensitivities with respect to the 
design parameters. 

The primary goal of this report is to describe results (Case 1 - Case 4) of 
simulations performed on different NexTech cell configurations. Also, development of 
fluid-structure interaction capability to analyze thermal stress distribution (Case 4, Case 5) 
and stress sensitivities (Case 5) in different components of SOFC is described. The results 
presented in Case 1 deals with the FlexCell, which contains honeycomb electrolyte. The 
fuel cell code required several modifications to analyze this unconventional geometry. The 
cell used in the second case (Case 2) can be termed as a simple fluid cell. Cell components 
included in this simulation are fuel manifold, current collector and interconnect. Primary 
purpose of this analysis was to compare results obtained using the SimCenter code with the 
numerical results obtained by the NexTech personnel. The third geometry (Case 3) 
obtained from NexTech was a 28mm cell that included all relevant components. However, 
while the grid was being generated on this cell, it was decided through mutual discussion 
to shift focus to the more important G13 cell. A thorough analysis on G13 cell (Case 4) has 
been performed that also includes stress analysis. Finally, results obtained on a simplified 
SOFC cell have been shown in Case 5. This geometry was utilized to develop initial fluid- 
structure interaction capability to perform stress analysis. This case also demonstrates the 
capability to compute sensitivity derivatives of structural cost functions. 

2. Solution Methodology 

2.1 Governing Equations (Fuel Cell) 

The three-dimensional SOFC model [30,31] utilized in this report solves the multi¬ 
species Navier-Stokes equations along with an electric potential equation that governs the 
distribution of electric potential and current density in the field. The three-dimensional 
model accounts for all components of the SOFC, including the anode, cathode, electrolyte, 
current collectors, seals, interconnects, and the fuel and air channels. Note that the model 
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is not limited to any particular type of SOFC, i.e. planar as well as tubular type SOFC can 
be simulated using this model. 

The governing equations for mass, momentum and energy conservation are solved 
simultaneously with the equation governing the electric potential in the numerical model. 
The system of equations utilized in the model is given by equations (1) - (6), which 
represent the conservation statements for the species concentrations, momentum (x, y and 
z), energy and current, respectively. 

Equations (1) - (5) are modified Navier-Stokes equations valid for both porous and 
fluid regions. Detailed discussion on flux formulation for these equations can be found in 
previous work [28-31]. Equation (6) represents the electric potential equation. As solid 
regions are considered zero-velocity regions, only energy (interconnect, electrolyte, seals) 
and electric potential equations (interconnect, electrolyte) are solved inside them. 
Electric/ionic conductivity, a, in equation (6) is a strong function of the temperature. 
Expressions describing the relationships between the electric/ionic resistivity (reciprocal of 
conductivity) and the temperature for various components of SOFC are presented in Table 
1 [44,45] along with thermal conductivities and other material properties of different 
components of the SOFC. 


d(sp.) ~ — 

+ V»(epT )+V»(J)-S 
at 


(1) 

d(epu) - dP e 2 up 

-f +V.(£pwC)= e +V*(£T ) 
at ox B 

(2) 

d(epv) _ / - dP _ . . £ : vp 

3 , +V-(£pvF)= £^,+V.(£T,) b 

(3) 

d(epw) - dP £ 2 H'P 

; ' + V.(£ pwV)= £ + V«(£T ) 

at az B 

(4) 

v<£(£ + P)V)+V'(Y7h.) = V'( £ prV) V-q*+V</>.(oV</>) 

dt ns 

(5) 

V.(oVp) = 0 

(6) 


As presented, equation (6) is an elliptic equation contrary to the rest of the 
governing equations, which are hyperbolic-parabolic equations. Equation (6) is solved in 
the entire domain except for the fuel and air channels, which are pure fluid regions. 

Several transport processes take place at the anode-electrolyte and the cathode- 
electrolyte interfaces that strongly affect the overall behavior of the SOFC. The 
electrochemical reactions taking place at the cathode-electrolyte and anode-electrolyte 
interfaces can be described by equations (7) and (8), respectively. 
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Table 1. Material properties of various components of SOFC [44, 45] 

Electric resistivity of anode (Qjw) 

2.98 x 10" 5 exp(-1392 / T) 

Electric resistivity of cathode 

8.11 x 10" 5 exp(600 / T) 

Electric resistivity of interconnect (flm) 

6.41 xlO' 8 

Ionic resistivity of electrolyte (Qm) 

2.94 x 1(T S exp(10350 / T) 

Thermal conductivity of anode (fV ) 

40.0 

Thermal conductivity of anode current collector W wT'AT -1 ) 

11.0 

Thermal conductivity of cathode [w m~' 

10.0 

Thermal conductivity of cathode current collector ( Wm~ x K~ x \ 

11.0 

Thermal conductivity of interconnect (w 

25.0 

Thermal conductivity of electrolyte (w 

2.0 

Thermal conductivity of seals [W m~ K ~ ] ) 

1.5 

Porosity of anode 

0.6 

Porosity of cathode 

0.6 

Tortuosity of anode 

6.0 

Tortuosity of cathode 

6.0 

Porosity of anode current collector 

0.9 

Porosity of cathode collector 

0.9 

Tortuosity of anode collector 

1.5 

Tortuosity of cathode collector 

1.5 


The effect of the aforementioned electrochemical reactions is modeled by applying 
mass flux conditions at the cathode-electrolyte (equation (9)) and anode-electrolyte 
(equations (10)-(11)) interfaces using Faraday’s law. 


J =- M, 

°2 4 f 0 


J„ = - M 

n 2 2F h 


(9) 

( 10 ) 
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( 11 ) 



In the above equations, i is the local current density and F is Faraday’s constant. 

A negative sign implies that the flux is leaving the interface. 

To account for the heat generated due to electrochemistry, heat flux proportional to 
the entropy change associated with the electrochemical reaction is applied at the anode- 
electrolyte and cathode-electrolyte interfaces. This heat flux is proportional to the molar 

formation rate, (i / n e Fj , where « is the number of electrons participating in the 

electrochemical reaction. 

In addition to the electrochemical reactions, two chemical reactions, methane 
reforming (12.1) and water gas shift (12.2) reactions have been included inside the anode 
for Case 5. For other SOFC cases (Case 1, Case 4), the fuel mixture contains only two 
species, hydrogen and steam and thus, only electrochemical reactions described by 
equations (7) — (8) arc included in the model. 


CH a + H^O <=± CO + 3 H, 

4 2 kb 2 


( 12 . 1 ) 


CO + r ± — > CO. + H. 

2 kb. 2 2 


( 12 . 2 ) 


Reaction rates for various species have been computed using the following 
equations. 


Rate r — kf r p CH4 P H2 o kb r p co p H1 
Rate s = kf,p co p H20 ~ kb s p CQ2 p H2 


(13.1) 

(13.2) 


Subscripts “r ” and “ j ” stand for reforming and shift reactions, respectively. 
Reaction rate constants, kf and kh , are computed using the methodology outlined in 
reference [29]. 

The voltage output of the SOFC strongly depends on several irreversibilities or 
losses encountered in the flowfield including activation polarization, concentration 
polarization and Ohmic polarization. Noren and Hoffman [46] have provided extensive 
discussion on accurately modeling the activation polarization. The SOFC model used in 
this work employs the Butler-Volmer equation to compute activation polarization [46]. 

The Butler-Volmer equation can be written as, 
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The activation polarization is denoted by i] icl . 

a is the charge transfer coefficient and assumed to be 0.5 in the current work. 
n e represents the number of electrons involved in the electrochemical reaction, which is 2 
(equations (7) and (8)) in the current simulation. 

i Q is the exchange current density and is computed using equations (15) and (16) for the 
anode and cathode [47], respectively. 


^0 ,a **a 

^0 ,c C 


D V 


P , , 

ref 7 


■ H,0 


P. , 

ref 7 


exp 


act,a 

R T 


P \ 025 

o. 


P , . 

ref 7 


exp 


act,c 

~RT , 

« 7 


(15) 

(16) 


Various constants in the above equations are given in Table 2 [47]. Once the 
values of a and n e arc inserted in equation (14), the activation polarization can be 
computed using the following expression. 


R T) 

u 

sinh -1 

( . \ 

/ 

F J 




(17) 


Table 2. Constants used to compute activation polarization [47] 

a 

0.5 

n 

e 

2 

ZM m- 2 ) 

5.5 xlO 8 

C r ( A m ~ 2 ) 

7.0 xlO 8 

E^Jhnor') 

1.0 xlO 8 

E oac (J kmol-') 

1.2 xlO 8 

PAN m- 2 ) 

101325 


Ohmic polarization is a direct consequence of the resistance offered to the flow of 
electrons/ions inside various components of the SOFC. Voltage drop due to Ohmic 
resistance is directly proportional to the current and the resistance. The effect of Ohmic 
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polarization on the voltage loss is directly included in the potential equation, equation (6), 
through the electric conductivity, <7 , which is the reciprocal of the electric resistivity. 

Concentration polarization is caused by reductions in the concentrations of the 
reacting species at the interface between the electrodes and the electrolyte. The effect of 
the reduction in concentrations can be seen from the well-known Nemst potential equation, 
given by equation (18). Also, exchange current densities at the anode-electrolyte interface 
and the cathode-electrolyte interface, represented by equations (15) and (16), respectively, 
are strongly affected by the concentration polarization. 


Equation (18) computes the electromotive force (EMF) or electric potential under 
reversible conditions, i.e. in the absence of activation, Ohmic or any other losses. 


RT 

EMF = EMF 0 +-In 

2 F 


( U pO-5 ^ 


V 


n,o 


- P 
.where P -—— 

1 p f 

ref 


(18) 


The electromotive force at standard pressure, EMF 0 , is computed using polynomial 
thermodynamic relationship between Gibbs free energy and temperature of different 
species participating in the electrochemical reactions. The value of P re) is taken as one 
atmosphere in the above equation. 


The electrochemical reaction reduces the concentration of the reactants and 
increases the concentration of the products at the electrode-electrolyte interface. Thus, the 
partial pressures of the reactants and products are affected in the same manner. This will 
reduce the value of the second term on the right-hand side of the equation (18) thereby 
affecting the EMF of the cell adversely. Concentration polarization strongly depends on 
the material properties of the electrodes that are responsible for the transport (diffusion and 
convection) of the reactants and products, to and from the electrode-electrolyte interface. 


2.2 Governing Equations and Computational Methodology (Structure) 

2.2.1 Governing Equations (Structure) 

As mentioned earlier, capability to perform thermo-mechanical stress analysis has 
been developed to analyze stress inside different components of the SOFC. Thermo¬ 
mechanical stress analysis can be performed either fully coupled, or one-way coupled in 
which only one disciplinary response affects the other. Regardless of the formulation, the 
governing equation may be expressed as 
d 2 u 


o +b. 


!/.' J 



(19) 


where, <7 represents the stress tensor, b the body force terms per unit mass, p the mass 
density, and u the displacement field. To cast the above equations in terms of 
displacements, the relationships between strain-displacement and stress-strain must be 
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assumed. For small strains and displacements (i.e., geometric linearity), the strains are 
related to the deformation gradients as 

£ <J = \{ U JJ + U u) ( 20) 

Additionally, under the assumption of linear elasticity (i.e. material linearity), the 
stress may be related to strain as 


a.. = C....(s.. -a-AT-8..) 

IJ ijU\ V kl I 


( 21 ) 


where, a is the coefficient of thermal expansion (which may in general be a function of 
temperature), AT the temperature difference from reference, and 8 the Kronecker delta 
symbol. Here, the first term relates the stress to mechanical strain, and the second to the 
thermal strain. Assuming isotropic material behavior, the constitutive (elasticity) tensor 
may be conveniently written as 


C 


ijkl 




Ev 

(l + v)(l-2v) 


8 . 8 ., 

ij kl 


( 22 ) 


where, in general, the modulus of elasticity, E, and Poisson’s ratio, V, may be functions of 
temperature. Currently, it is assumed that the mechanical response does not alter the 
temperature distribution, and therefore, a one-way coupling is utilized. Furthermore, a 
steady-state temperature field is applied to the structure and, hence, inertia may be 
neglected in the problem formulation. 


2.2.2 Computational Methodology (Structure) 

The thermo-elastic structural analysis is performed using a standard displacement- 
based Galerkin formulation. With introduction of the stress-strain, and strain-displacement, 
relations into the equations of equilibrium, the Navier-displacement equations may be 
written (neglecting inertia) as 


<V.,„ +i /=° < 23 ) 

Integrating these equations over the volume of an element, using a standard 
Galerkin formulation, and assembling the element equations yields an algebraic system to 
be solved for the unknown nodal displacement vector {d} as 


where, the global stiffness matrix and load vector are 

M = XJW[ C ]I>!K (25) 

{^} = I|[v]’{(.}rfF+XJ[v] r {(}j5+Xj[A'] r [c]{a.Ar}rfF (26) 


Here, ("a! represents the element shape function matrix, which relates the 
displacement field to the element nodal displacements, and is the so-called strain- 
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displacement matrix, which relates the element strains to the nodal displacements. The first 
term of the load vector represents body forces acting over the volume of the domain, the 
second term represents forces due to traction or stress (t. - 0.7] , where 7). are the 

components of the unit outward pointing normal) acting on the boundaries of the domain, 
whereas the last term are the loads resulting from the thermal strains. 

The solution to the resulting system may be accomplished with either an iterative or 
direct method. Iterative solution algorithms can be beneficial for very large systems which 
may prohibit direct solution methods, or for systems which do not have symmetric 
damping matrices. The direct method uses a sparse Cholesky decomposition based on 
skyline storage scheme. A preconditioned GMRES [48] is utilized to solve the system 
iteratively. 

Once the nodal displacements have been determined, the strain and stress fields 
may be computed. In multidimensional problems, the individual components of the stress 
tensor do not provide adequate information in order to ascertain the proximity to failure or 
yielding of the material. In these regards, different stress measures are typically used. For 
brittle materials, the maximum in plane principal stress is typically used. The principal 
planes are those in which the shear stress vanishes and, therefore, are the planes that have 
maximum normal stress. Since brittle materials tend to fail due to normal stress, 
appropriate failure criteria are usually based on maximum principal stress. For ductile 
materials, which tend to fail in shear, typically either the von Mises or the Tresca yield 
criteria are used. The von Mises yield criterion uses the assumption that the onset of yield 
is based on the second deviatoric stress invariant. The Tresca yield surface is 
circumscribed by the von Mises yield surface, representing a more conservative criterion 
for prediction of plastic yielding. 

2.3 Boundary Conditions 

Boundary conditions utilized in the numerical solver are described for a sample 
geometry shown in Figure 1. The geometry includes all relevant SOFC components 
including air/fuel channels, anode, cathode, electrolyte and interconnects. No-slip, 
adiabatic wall boundary conditions are applied at the interfaces between the electrodes and 
the interconnect, as well as the side, top and bottom walls. Solid regions (interconnect and 
electrolyte) are considered zero-velocity regions and the only variables computed inside 
these regions are temperature and electric potential. A fixed potential (0 = 0) boundary 
condition is applied at the bottom wall, whereas the top wall is treated by specifying 
average current density (i = i Ued )• 

Inflow boundary conditions with specified mass flow rate and species mole 
fractions are applied at both fuel and air channel inlets. Specified back pressure outflow 
conditions are applied at both air and fuel channel outlets. 
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Figure 1. Boundary and interface conditions (I - Interconnect, C - Cathode, AC - Air 
Channel, E - Electrolyte, FC - Fuel Channel, A - Anode) 


2.4 Sensitivity Analysis - Direct Differentiation 

As mentioned earlier, direct differentiation method is utilized in both the fuel cell 
and structures codes to compute sensitivity derivatives. Derivation of this method using the 
chain rule is shown in equations (27) - (30). Q represents the vector of solution variables 
where each element of the vector is representative of one or more physical variables 
located at each mesh point, X ■ R represents the vector of discrete residuals at each mesh 
point, / is the cost function and p is a vector of design variables. 

dp dp dQ dp dx dp 

Now, R(Q(P),xUB),P)-0 

M + ^_dQ + Mdx _ Q 

dp ~ dp dQ dp dx dp ~ 


(27) 

(28) 
(29) 
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As seen, computation of dQ / dp is an essential component of this method, which 
requires the solution of a linear system of equations for each design variable. This 
requirement makes direct differentiation methods computationally expensive for problems 
with many design variables. Discrete adjoint method is more favorable for the problems 
involving large number of design variables. However, as derivatives of dependent 
variables with respect to the design variables are computed at each node in the flowfield in 
direct differentiation, this method is particularly useful when there are many flowfield 
constraints. 


2.5 Solution Procedure 

Flowfield variables are computed using an unstructured, implicit, finite-volume 
solver. The solver is vertex centered and the discrete residual at each node is computed by 
integrating the governing equations, (1) - (6) over a median dual control volume. Because 
a steady-state solution is the primary goal of the current work, time accuracy of the 
solution is sacrificed by allowing local time-stepping to accelerate convergence. 

To reduce computer time, the solution is obtained using multiple processors 
utilizing the message passing interface (MPI) [49] and necessary grid decomposition is 
achieved using METIS [50]. Original grids are generated using the commercial software 
Pointwise [51]. 

An implicit Euler scheme is used to solve the non-linear system as given by 
equations (1) - (6). A flux-difference splitting scheme based on the ROE scheme [52,53] 
for a multi-component mixture is derived to model the convective fluxes. A central- 
difference formulation is used to compute all the second-order derivative terms. Linear 
systems encountered in both the flowfield and sensitivity solvers are solved using the 
GMRES [48] method. 

3. Results and Discussion 

3.1 Analysis 

3.1. Case 1 (FlexCell) 

The first NexTech configuration analyzed in this study was a FlexCell. Dimensions 
of the various components of the FlexCell were obtained from NexTech along with two 
separate CAD files of the fuel manifold and honeycomb electrolyte. Significant efforts 
were required to integrate these CAD models with the anode and cathode layers to obtain a 
complete CAD model for the FlexCell. After the development of the CAD model, an 
incremental approach was adopted to perform numerical simulations by progressively 
adding more complexity. Initial simulations included the fuel manifold geometry without 
accounting for the electrochemistry. The aim of these simulations was to observe the flow 
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distribution among different channels. Because the honeycomb electrolyte of the FlexCell 
poses several numerical challenges, several simulations included a flat electrolyte and 
channels that were also shortened and reduced in number to lessen the number of grid 
elements, which are of the order of few millions. Finally, simulations were conducted that 
included honeycomb electrolyte. However, to reduce the number of elements the size of 
the geometry was reduced by one fourth and the plenum regions for both air and fuel 
manifolds were removed for simplification. One of the major challenges of this particular 
simulation was to generate a good quality grid, especially for the honeycomb regions. In 
our previous work, and most other numerical simulations reported by researchers, the 
anode-electrolyte interface shape has been planar. Thus, another challenging task was to 
handle the irregular shape of the interface between the anode and the electrolyte. 

Figure 2(a) shows the flowfield geometry of the FlexCell. The manifold shapes are 
identical for both the fuel and air sides portions of the cell. Circular cross-sections shown 
in the figure indicate flow inlet and outlet ports. Because the initial simulations were 
conducted to analyze the flowfield in different channels of the manifold, 
electrochemistry/chemistry was not included. The velocity distribution on a plane passing 
through the fuel manifold is shown in Figure 2(b). 



Figure 3(a) shows different layers of the FlexCell. As seen, the anode-electrolyte 
interface is not flat, which adds significant complexity in the numerical simulations. Also, 
generating grid for honeycomb regions is extremely difficult as shown in Figure 3(b). As 
seen, each hexagon is broken into several regions and sizes/shapes of various sub-regions 
are different for each hexagon. 
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Figure 3(a). FlexCell layers 


Figure 3(b). Grid methodology 


Figure 3. Different layers of the FlexCell and grid methodology 


Anode Layer 
Barrier Layer 


Electrolyte Support Layer 


Electrolyte Membrane Layet 
Barrier Layer 
Cathode Layer 


Figure 2(i). FlexCell Layera 


Figures 4(a) and 4(b) show results obtained on a truncated FlexCell with all 
relevant transport processes and electrochemistry included in the simulation. Mole 
fractions of different species included in this simulation along with operating conditions 
are described in Table 3. To control the number of grid elements, the geometry of the 
FlexCell has been reduced by one-fourth and the plenum regions for both fuel and air 
manifolds are removed. 


Table 3. Mole fractions and thermodynamic conditions 

** 

X„,o 


** 

T{K) 

P(N/m 2 ) 

0.95 

0.05 

0.21 

0.79 

1073 K 

101325 


Figure 4(a) shows temperature contours in a plane passing through the honeycomb 
region. Because the electrochemical reaction associated with an SOFC is exothermic, a 
gradual increase in temperature is evident in the flow direction. The layer shown in the 
figure is made up of two different materials, namely anode and electrolyte. Both these 
materials have different thermal conductivities, which is responsible for the rough nature 
of the temperature contours shown in Figure 4(a). 
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Figure 4(b) shows steam concentration in a plane passing through the honeycomb 
region. Arrows in the figure indicate flow direction. Blue color in the figure indicates the 
electrolyte regions, which do not contain any chemical species. During the SOFC 
operation, steam is generated due to the electrochemical reaction. Thus, concentration of 
steam can be seen increasing in the figure. However, increase in steam concentration is not 
uniform in different hexagon regions. Various transport processes involved in SOFC 
operations and the location of each honeycomb with respect to the fuel channels are 
responsible for such behavior. 




Figure 4(a). Temperature contours (K) 


Figure 4(b). Steam concentration 


Figure 4, FlexCcll: Results of numerical simulations 


3.2 Case 2 (Simple Fluid Cell) 

The reason behind specifying this cell as a “simple fluid cell” is that the analysis of 
this cell only includes fluid flow in different components of the cell. The primary purpose 
of this study was to compare numerical results obtained using the SimCcntcr code with the 
simulation results obtained by NexTech personnel. Different components of simple fluid 
cell are shown in Figure 5, which includes interconnect, flow manifold and current 
collector. 



t 


Interconnect 

Flow manifold 

Current collector 

Figure 5. CAD model of simple fluic 

cell 
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The current collector utilized in this model is made up of highly porous material 
(porosity = 0.9). As can be seen in Figure 5, fluid is entering vertically at the inlet port at 
800°C and atmospheric pressure and flows through the current collector before leaving 
through the outlet port. The fluid containing the mixture of hydrogen and nitrogen with 
equal mole fractions is utilized in this simulation. 




Figure 6(a). Pressure distribution 


Figure 6(b). Contours of velocity magnitude 


_ Figure 6(c). Velocity vectors _ 

Figure 6. Simple fluid cell: Results of numerical simulations 


Velocity Mug nit ud 
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Figure 6(a) and (b) show pressure and velocity magnitude contours plotted on a 
plane extracted through the flowfield in streamwise direction. As expected, pressure is 
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reduced as fluid moves from the inlet towards the outlet port. NexTech personnel 
successfully compared the pressure contours plotted in figure 6(a) with their in-house 
numerical results both quantitatively and qualitatively. Figure 6(b) shows contours of 
velocity magnitude. The figure shows some interesting characteristics including, strong 
acceleration as flow enters from the inlet port to the main flowfield through two narrow 
passages and presence of regions with extremely low velocity magnitude located at four 
comers of the flowfield. Results shown in Figure 6(b) also compared favorably with the 
NexTech results. To further demonstrate the flow patterns, velocity vectors are plotted on a 
streamwise plane in Figure 6(c). As seen, vectors show uniform pattern throughout the 
domain except in the comers. 

3.3 Case 3 (28 mm cell) 

The third CAD model received from NexTech is shown in Figure 7. The 
computational model contains all SOFC components including flow channels, electrodes, 
electrolyte, interconnects, current collectors and seals. Figure 8(a) and 8(b) show shapes of 
anode and cathode current collectors, respectively. As seen, geometry of this cell is quite 
complicated mainly due to the shape of the current collectors and the manner in which 
current collectors are connected with other components of the cell. While model was still 
in the meshing phase, the decision to move to the more important G13 cell was made 
during discussion with NexTech. Even though analysis was not performed on this cell, 
significant efforts were made to generate mesh for this complicated computational model. 



Figure 7, CAD model of 28 mm NexTech cell 
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3.4 Case 4 (G13 cell) 

As mentioned earlier, analysis of G13 cell was prioritized over 28 mm cell after 
discussion with NexTech personnel. Different components of G13 cell are shown in 
Figures 9 and 10. Fuel flowfield, made up of fuel manifold and anode current collector is 
shown in Figure 9(a). The design of the fuel flowfield is inspired by the simple fluid cell 
described in Case 2. Figure 9(b) shows air flowfield comprising of air manifold and 
cathode current collector. In the initial CAD model received from NexTech, shapes of the 
anode and cathode currents collectors were not exactly the same. However, it was later 
decided to use the same shape (as the anode current collector) for both current collectors to 
simplify the model. Electrodes and electrolyte are sandwiched between the flowfields 
shown in Figure 9. To avoid mixing of fuel and air, seals are utilized on both sides (Figure 
10(a)). All these components are enclosed on top and bottom by the interconnects shown in 
Figure 10(b). 



The inlet and outlet ports for both air and fuel manifolds are indicated in Figure 9. 
As air and fuel flow in the opposite directions, the flow arrangement is counter-flow. Fuel 
mixture and air enter through appropriate inlet ports and diffuse inside current collectors 
and electrodes. Oxygen atoms combine with the electrons inside an extremely thin layer 
near the cathode-electrolyte interface and are converted into oxygen ions 
(O.50 2 + 2e~ —> 0 2 ~). These oxygen ions migrate through the solid electrolyte to reach the 
anode-electrolyte interface and combine with hydrogen to generate steam and release 
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electrons (//. + O 2 ' -* H 2 0 + 2e ). Both these reactions take place inside extremely thin 

layers near the cathode-electrolyte and the anode-electrolyte interfaces, respectively. 
Because modeling the details of the interface region is impractical due to the small size, the 
cumulative effect of the electrochemical reactions is modeled as a jump in the electric 
potential, which is determined using Nemst potential (equation (18)) and activation 
polarization (equation (17)). 



The fuel mixture is assumed to contain hydrogen and steam only, i.e. chemistry' is 
not present in this case. Air is modeled as a mixture of oxygen and nitrogen. Species mole 
fractions of the fuel mixture and air entering the flowficlds arc given in Table 4 along with 
the operating pressure and temperature of both fuel and air. 


Table 4. Mole fractions and thermodynamic conditions 


*h 2 o 


** 

T(K ) 

P(N /m 1 ) 

0.95 

0.05 

0.21 

0.79 

1073 K 

101325 


Figures 11-13 show temperature contours plotted on planes extracted through 
different components of the G13 cell in streamwise direction. Temperature is continuous in 
the entire domain unlike potential or concentration fields. As mentioned earlier, various 
modes of energy transfer includes convcction/diffusion/conduction of energy, heat 
generated due to viscous stresses and most importantly, heat generated due to entropy 
changes occurring because of the electrochemical reaction at the anode-electrolyte 
interface. Also, Ohmic heating contributes to the rise in temperature. 

The extent of temperature rise depends strongly on the mass flow rates of fuel and 
air. Convective cooling increases with higher flow rates of either fuel or air and thus, 
contributes towards reduction in average temperature rise. The opposite is true for the 
lower mass flow rates. Current density is another major contributor in temperature rise as 
Ohmic heating is linearly proportional to the current density. Also, current density dictates 
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the extent of electrochemical activity at the electrode-electrolyte interface that is 
responsible for generating heat. 


Figures 11(a) and (b) show temperature contours plotted on planes passing through 
fuel and air manifolds, respectively. Flow directions of both fuel and air are also indicated. 
In figure 11(a), fuel enters the inlet port at 1073 K and there is a gradual rise in fuel 
temperature due to the heat generated by electrochemical reaction. As fuel moves through 
the flowfield, temperature starts reducing due to the cooling effect of the incoming air from 
the opposite direction. The reduction in fuel temperature continues untill it leaves through 
the outlet port. In figure 11(b), air enters through the inlet port at 1073 K and gradually 
heats up as it passes through the flowfield. The rise in the temperature of air is continuous 
throughout the domain unlike fuel until it reaches the outlet port, when it begins to cool 
down. Overall, temperature of fuel is lower than the temperature of air leaving from the 
respective outlets. 


Temperature (K) 



0 



Temperature (K) 

1145 
1140 
1135 
1130 
1125 
1120 
1115 
1110 
1105 
1100 
1095 
1090 
1085 
1080 
1075 


Figure 11(a). Fuel manifold 


Figure 11(b). Air manifold 


Figure 11. Temperature contours 


Figures 12(a) and (b) show temperature contours plotted on planes passing through 
the anode and cathode current collectors, respectively. A lot of similarities can be found 
between the characteristics of these contours and the contours shown in figure 11 (a) and 

(b) including spots indicating higher and lower temperature regions. Figures 13(a), (b) and 

(c) show temperature contours plotted on planes passing through the anode seal, cathode 
seal and electrolyte, respectively. Again, temperature distribution inside these components 
exhibits similarities with those shown in Figures 11 and 12. Similarities in temperature 
contours inside different components of the cell indicate that the temperature gradient 
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through the thickness of the cell is negligible due to the small thickness of the cell 
compared to other dimensions. 




Figures 14(a) and (b) show hydrogen and steam contours plotted on a plane passing 
through the fuel manifold in streamwise direction. Figure also indicates the flow direction 
of the fuel. In electrochemical reaction, hydrogen acts as a reactant while steam is a 
product. Thus, hydrogen concentration reduces as fuel flows from the inlet port towards 
the outlet port. Non-uniformity in reduction of hydrogen concentration is evident in Figure 
14(a), especially left and right sides of the domain show the lowest concentration. One of 
the reasons behind this behavior is that the fuel flow rate in the affected area is low 
compared to the rest of the domain. Further discussion regarding flow rate is presented in 
the later part of this report. On the other hand, steam concentration, which is a product of 
the electrochemical reaction, increases as fuel flows through the domain as shown in 
Figure 14(b). 
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Figure 14(a). Hydrogen concentration 
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Figure 14(b). Steam concentration 


Figure 14. Species concentration in fuel manifold 


Figures 15 and 16 show contours of hydrogen and steam concentration on planes 
extracted in streamwise direction through anode and anode current collector, respectively. 
Characteristics of concentration contours in figure 15 and 16 are similar to those discussed 
for figure 14 in previous paragraph. This suggests strong species diffusion through the 
thickness of the cell. One of the reasons for such behavior is small thickness of different 
components of the cell. 


rho_H2(kg/m3) 



rho_H20(kg/m3) 



Figure 15(a). Hydrogen concentration 


Figure 15(b). Steam concentration 


Figure 15. Species concentration in anode 
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Figure 16(a). Hydrogen concentration 

Figure 16(b). Steam concentration 

Figure 16. Species concentration in anode current collector 


Figures 17(a), (b) and (c) show contours of oxygen concentration on planes 
extracted through air manifold, cathode and cathode current collector, respectively. As 
oxygen is a reactant of the electrochemical reaction, its concentration reduces as air flows 
through the cell. As shown for the hydrogen and steam, oxygen contours plotted through 
different components in figure 17 exhibit similar characteristics indicating strong species 
diffusion through the thickness of the cell. Spots with the lowest oxygen concentration 
(top-left and top-right comers) have strong relationship with the air flow rate. Further 
discussion on this topic is covered in the later part of this report. 
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_ Figure 17(c). Cathode current collector _ 

Figure 17. Oxygen concentration inside different components of G13 


Figures 18(a) and (b) show contours of streamwise velocity (z-direction) plotted on 
planes passing through fuel and air manifolds, respectively. These are some interesting 
features present in both plots. Non-uniformity in fuel flow rate is evident in figure 18(a) as 
majority of fuel is passing through the middle region of the manifold. The reasons behind 
such behavior include the placement of the inlet port with respect to the main flowfield and 
the angles of the narrow passage through which fuel enters the flowfield. The starving 
areas are located at four comers of the flowfield in figure 18(a). Unlike fuel manifold that 
has one inlet and exit port, air manifold contains two inlet and exit ports. Thus, flow 
distribution inside air manifold is completely different that that of the fuel manifold. Due 
to the placement of the inlet ports and angle at which air inters the main flowfield, areas 
with high flow rate are located between the center and the sides of the domain in figure 
18(b). Again, the starving regions are located at the comers in figure 18(b). For both air 
and fuel manifolds, comers represent the most affected regions due to low flow rates of 
respective fluids, which may adversely affect the cell efficiency and lead to localized 
heating. 

Figures 19-23 show stress distribution inside different components of the G13 cell. 
All components are considered stress free at 800°C. The stress contours are generated by 
the in-house structures code using flowfield data (temperature) imported from the fuel cell 
code. SOFC is manufactured by connecting various components that possess different 
material properties. Thus, mismatch in coefficients of thermal expansion of these 
components contribute greatly towards stress generation in SOFC. Figures 19(a) and (b) 
show contours of mean principal stress (MPS) plotted on outer surface of the anode and 
cathode current collectors, respectively. As seen, stress distribution in both current 
collectors is uniform except at some outer edges, where current collectors are connected 
with other cell components. 
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MPS (MPa): *350 *256 *162 -68 26 120 



MPS (MPa): -050 -350 -250 -150 *50 50 150 


Figure 19(a). Anode current collector [ Figure 19(b). Cathode current collector 

_ Figure 19. Stress distribution _ 


Figure 20 shows MPS contours plotted on a plane extracted through the electrolyte 
in streamwise direction. Higher values of MPS are present near the active area, where 
electrochemical reaction takes place. Also, electrolyte is connected with anode and cathode 
on each side of this region. Thus, mismatch in coefficients of thermal expansion and 
presence of exothermic electrochemical reactions are responsible factors for high stress in 
this region. 










































Figure 20, Stress distribution inside electrolyte 
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Figures 21(a) and (b) show maximum tresca equivalent stress (MTES) plotted on 
planes passing through seals located on fuel and air sides, respectively. Overall, stress 
distribution is uniform in both figures. A close observation indicates isolated spots with 
high stress values located on the edges surrounding inlet and outlet ports of both fuel and 
air manifolds. 


Figures 22(a) and (b) show MTES contours plotted on planes passing through the 
interconnects located at the bottom (fuel side) and top (air side) of the cell, respectively. As 
clearly seen, regions with maximum stress are located near air outlet ports in both figures. 
Also, tiny spots with high stress values are present near fuel inlet port in both figures. 
Characteristics of stress contours in Figure 22 are similar to those for the seals shown in 
Figure 21. 



Figure 22(a). Interconnect (fuel side) 

Figure 22(b). Interconnect (air side) 

Figure 22. Stress distribution 


3.5 Case 5 (Simplified SOFC) 

The fluid-structure interaction capability demonstrated in the previous case was 
initially developed for the simplified SOFC described in this case. The development was 
performed while waiting for the arrival of new cluster and storage system to be used to run 
simulations on NexTech’s G13 cell. The geometry of the cell utilized in this case is shown 
in Figure 23. As seen, computational model includes all relevant SOFC components 
including fuel/air manifolds, electrodes/electrolyte (PEN), seals and interconnects. The 
number of channels for both air and fuel manifolds is chosen to be twelve to make the 
geometry mimic realistic planar SOFC. Various geometric dimensions of the cell are listed 
in Table 5. 
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Table 5. Geometrical dimensions of simplified SOFC 

Length 

43.0 mm 

Width 

39.0 mm 

Height 

0.23 mm 

Anode thickness 

0.1 mm 

Cathode thickness 

0.1 mm 

Seal thickness 

0.1 mm 

Electrolyte thickness 

0.1 mm 

Interconnect thickness 

0.5 mm 

Channel thickness 

0.5 mm 


No-slip, adiabatic wall boundary conditions are applied at the top wall, bottom wall 
and side walls of the computational geometry shown in Figure 23. As mentioned 
previously, fixed potential (0 = 0) boundary condition is applied at the bottom wall, while 
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the top wall is treated by specifying average current density (/ = i applied ). Inflow boundary 

conditions with specified mass flow rate and species mole fractions are applied at both fuel 
and air channel inlets. The temperature of the air and fuel mixture entering from their 
respective inlet ports is 1073 K. Also, both fluids are operating at atmospheric pressure. 
Specified back pressure outflow conditions are applied at both air and fuel outlet ports. 
Initial species mole fractions and thermodynamic conditions utilized in this simulation are 
given in Table 6. As seen in Table 6, partially reformed fuel has been utilized and thus, 
methane reforming and water gas shift reactions have been considered inside the anode. 
Current density of 5500 Am' 2 is applied at the top wall of the computational geometry 
shown in Figure 23. 


Table 6. Mole fractions and thermodynamic conditions 

*co 


* C02 

** 

X CH , 

** 

** 

T(K) 

P(N m~ 2 ) 

0.029 

0.493 

0.044 

0.263 

0.171 

0.198 

0.802 

1073 K 

101325 


Two different configurations co-flow and counter-flow, are analyzed in this case. 
Figures 24(a) and (b) show temperature contours plotted over outer surface of the cell for 
co-flow and counter-flow configurations, respectively. Figures also show flow directions 
of both air and fuel. As expected, in co-flow case, there is a gradual rise in temperature as 
both fuel and air move through the flow domain. Heat generated due to the electrochemical 
reaction is the main factor affecting the increase in temperature. In counter-flow case, 
regions showing maximum temperature are present in the middle of the computational 
domain. Also, maximum temperature found in the co-flow case is higher than the same 
found in the counter-flow case. 




Temperature (K): 1060 HOO 1120 n*0 nfiO n80 1200 1220 1240 

Temperature (K); 1OB0 1100 1120 1>40 1160 11B0 130Q 



Figure 24(a). Co-flow configuration 

Figure 24(b)* Counter-flow configuration 

Figure 24. Temperature contours on outer surface on the cell 


Figure 25 shows polarization curves plotted for both co-flow and counter-flow 
cases operating under the same conditions described in Table 6. As expected, cell voltage 
reduces with increase in current density due to the effects of several irreversibilities present 
inside the cell. Both cases exhibit similar performance for low current densities. However, 
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as current density increases, co-flow configuration performs better than the counter-flow 
configuration. 



Figures 26 (a) - (f) show mole fractions of different species plotted on planes 
passing through fuel and air manifolds for the co-flow case. As mentioned earlier, tw r o 
chemical reactions namely, methane reforming and water-gas shift reactions are considered 
inside the anode electrode. Also, electrochemical reaction, which is responsible for the 
production of steam and consumption of hydrogen and oxygen, affects species distribution 
in the flowfield. In figure 26(a), there is an overall reduction in hydrogen concentration as 
it moves through the flowfield. A region located near bottom right comer of the plane 
shows rise in hydrogen mole fraction. This behavior is caused by hydrogen production due 
to methane reforming reaction. In figure 26(b), gradual rise in steam concentration due to 
electrochemical reaction is evident as fuel moves from the inlet to the outlet of the 
manifold. Methane mole fraction is plotted in figure 26(c). As methane reforming is a fast 
reaction, most of the methane can be seen consumed in the first half of the flowfield. 

Figure 26(d) shows oxygen mole fraction plotted on a plane extracted from the air 
manifold. As oxygen is a reactant of the electrochemical reaction, there is gradual 
reduction in its mole fraction as air moves through the flowfield. Contours of carbon 
monoxide (CO) mole fraction plotted in figure 26(e) exhibit non-uniformity in the 
flowfield. As CO acts as a reactant in shift reaction and as a product in reforming reaction, 
their combined effect produces non-uniformity in the contours shown in figure 26(e). 
Finally, contours of carbon dioxide mole fraction are plotted in figure 26(f). As the only 
reaction involving carbon dioxide (as a product) is a shift reaction, gradual rise in carbon 
dioxide concentration is evident in figure 26(f). 
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(d). Oxygen 


(e). Carbon monoxide |(f). Carbon dioxide 


Figure 26. Mole fraction of different species in fuel and air manifolds (Co-flow) 


Figures 27(a) - (f) show stress contours plotted on planes extracted through 
different solid and porous components of the cell for co-flow configuration. Figures 27(a) 
and (b) show contours of MTES plotted on streamwise planes passing through seals 
located on both fuel and air sides, respectively. As seen, maximum stress is present near 
fuel inlet port in both figures. Figure 27(c) - (e) show contours of MPS plotted on planes 
passing through the anode, electrolyte and cathode, respectively. Some characteristics of 
these contours such as location of the maximum stress region are similar in all three plots. 
This region is located near fuel inlet port in the computational domain. Figure 27(f) shows 
contours of MTES plotted on a vertical plane passing through the interconnect and inlet 
ports. As expected, maximum stress is found near fuel inlet port. Regions near air inlet port 
also indicate high values of MTES. Overall, maximum stress values are found in a plane 
extracted through the anode electrode in figure 27(c). 
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(a). Anode seal 



(c). Anode 



(d). Electrolyte 


b). Cathode seal 



(e). Cathode 


MTES(MPa): 30 50 70 90 110 130 


(f). Interconnect 


Figure 27. Stress distribution inside different components for co-flow configuration 


Figures 28(a) - (f) show stress contours plotted on planes extracted through 
different components of the cell for counter-flow configuration. Even though air is flowing 
in the opposite direction in this case compared to the previous case (figure 27), stress 
contours in different components show similar characteristics in both cases. Regions with 
maximum stress are located near fuel inlet port in different cell components. Overall, stress 
values for the counter-flow case are smaller than the co-flow case. 

As mentioned earlier, both fuel cell and structures code are capable of performing 
sensitivity analysis. To demonstrate this capability, stress sensitivity contours are plotted 
on planes extracted through different components of the cell in figures 29(a) - (e). The co¬ 
flow configuration is utilized in figure 29. Design variable in this study is cathode porosity. 
The method utilized to compute sensitivity derivatives is direct differentiation in both fuel 
cell and structures code. Structures code requires values of flowfield variables and 
sensitivities of flowfield variables from the fuel cell code to compute stress sensitivities. 
The characteristics of stress sensitivity contours in figure 29 are similar to those shown for 
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the stress contours in figure 27; especially regions with highest sensitivity values are 
located near fuel inlet port for all components. 



MTES(MPa): 10 30 50 70 90 110 130 



(f). Interconnect 


Figure 28. Stress distribution inside different components for counter-flow configuration 
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(c). Anode 


(d). Electrolyte 


(e). Cathode 


Figure 29. Stress sensitivities with respect to the cathode porosity inside different 
components for co-flow configuration 


4. Conclusions 

A three-dimensional, parallel, unstructured solver has been developed to model 
complicated transport phenomena present inside all components (channels, electrodes, 
electrolyte, seals, current collectors and interconnects) of different cell configurations 
provided by NexTech Materials. A capability to perform thermo-mechanical analysis of 
different components of SOFCs has been developed by coupling the fuel cell code with the 
structures code. Results obtained using this capability indicate that the main factors 
affecting the stress distribution are temperature gradients and mismatch of coefficients of 
thermal expansion between different components of the cell. Also, capability to perform 
thermo-mechanical sensitivity analysis has been developed using direct differentiation 
method. Using this capability, sensitivity derivatives of stress has been computed for a 
simplified SOFC geometry. 
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Nomenclature 


Symbols 

Name 

Unit 

B 

Permeability 

m 2 

E, 

Total energy 

J m' 3 

f 

Cost function 

cost function depended 

H 

Enthalpy 

J kg' 1 

i 

Current density 

Am’ 2 

*0 

Exchange current density 

Am' 2 
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J 

Mass flux vector 

kg m' 2 s" 1 

• 

Mass flow rate 

kg s' 1 

m 


kg kmol’ 1 

M 

Molecular weight 

ns 

Number of species 

- 

P 

Pressure 

Nm’ 2 

Q 

Solution vector 

solution variable dependent 

q 

Heat flux 

J m' 2 s’ 1 

T 

Temperature 

K 

u 

x-velocity component 

m s' 1 

V 

y-velocity component 

m s' 1 

w 

z-velocity component 

m s" 1 

X, 

Mole fraction of i th species 

- 

Greek Symbols 



fi 

Design variable vector 

design variable dependent 

X 

Grid vector 

m 

P 

Mass concentration 

kgm' 3 

P 

Molecular viscosity 

kg m s’ 1 

<P 

Electric potential 

volt 

n 

Activation polarization 

volt 

£ 

Porosity 

- 

K 

Totruosity 

- 

V 

Control volume 

m 3 

T 

Viscous flux 

kg m" 1 s’ 2 

Constants 



F 

Faraday constant 

96484.56 A s mol’ 1 

K 

Universal gas Constant 

8314.4 J kmol' 1 K' 1 

Indices 



a 

Anode 
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c 


Cathode 


eff Effective 

i,j Chemical species 


References 

[1] Ferguson, J. R., Fiard, J. M., Herbin, R., “Three-Dimensional Numerical Simulations 
for Various Geometries of Solid Oxide Fuel Cells”, J. Power Sources, 58 (1996) 109-122. 

[2] Fiard. J. M., Herbin, R., “Comparison Between Finite Volume and Finite Element 
Methods for an Elliptic System Arising in Electrochemical Engineering”, Comp. Methods 
in App. Mech. and Eng., 115, (1994) 315-338. 

[3] Haberman, B. A., Young, J. B., “Three-Dimensional Simulation of Chemically 
Reacting Gas Flows in the Porous Support Structure of an Integrated-Planar Solid Oxide 
Fuel Cell”, Int. J. of Heat and Mass Transfer, 47, (2004) 3617-3629. 

[4] Khaleel, M. A., Lin, Z., Singh, P., Surdoval, W., Collin, D., "A finite element analysis 
modeling tool for solid oxide fuel cell development: coupled electrochemistry, thermal and 
flow analysis in MARC", J. of Power Sources, 130, (2004) 136-148. 

[5] Lehnert W., Meusinger J., Thom F., “Modelling of gas transport phenomena in SOFC 
anodes”, J. Power Sources, 87, (2000) 57-63. 

[6] Li P. W., Chyu. M. K., “Simulation of the Chemical/Elcctrochcmical Reaction and 
Heat/Mass Transfer for a Tubular SOFC Working in a Stack”, J. Power Sources, 124, 
(2003) 487-498. 

[7] Chan, S. H., Khor, K. A., Xia. Z. T., “A complete polarization model of a solid oxide 
fuel cell and its sensitivity to the change of cell component thickness”, J. Power Sources, 
93,(2001)130-140. 

[8] Campanari, S., Iora, P., J. Power Sources, “Definition and sensitivity analysis of a 
finite volume SOFC model for a tubular cell geometry”, 132 (2004) 113-126. 

[9] Grujicic, M., Chittajallu, K. M., “Design and optimization of polymer electrolyte 
membrane (PEM) fuel cells”, Applied Surface Sci., 227 (2004) 56-72. 

[10] Grujicic, M., and Chittajallu, K. M., “Optimization of the cathode geometry in 
polymer electrolyte membrane (PEM) fuel cells”, Chem. Eng. Sci., 59 (2004) 5883-5895. 


36 


[11] Anderson, W.K., Newman, J.C., Whitfield, D.L., Eric Nielsen, E.J., “Sensitivity 
Analysis for the Navier-Stokes Equations on Unstructured Meshes Using Complex 
Variables”, AIAA Journal,39, 1, (2001) 56-63. 

[12] Burden, R. L., Faires, D. J., “Numerical Analysis”, Sixth Edition, Brooks/Cole 
Publishing Company, 1997, 16-26. 

[13] Secanell, M., Djilali, N., Suleman, A., “Optimization of a planar self-breathing PEM 
fuel cell cathode”, AIAA 2006-6917, 11th AIAA/ISSMO Multidisciplinary Analysis and 
Optimization Conference, 6-8 Sept., 2006, Portsmouth, Virginia. 

[14] Huang, C. M., Shy, S. S., Lee, C. H., “On flow' uniformity in various interconnects 
and its influence on cell performance of planar SOFC”, J. Power Sources, 183, (2008) 205- 
213. 

[15] Lee, S., Jeong, IL, Ahn, B., Lim, T., Son, Y., “Parametric study of the channel design 
at the bipolar plate in PEMFC performances”, I. J. of Hydrogen Energy, 33, (2008) 5691- 
5696. 

[16] Anderson, W. K., Vcnkatakrishnan, V., “Aerodynamic Design Optimization on 
Unstructured Grids with a Continuous Adjoint Formulation”, Computers and Fluids, 28, 4- 
5,(1999) 443-480. 

[17] Anderson, W r . K., Bonhaus, D. L., “Airfoil Design on Unstructured Grids for 
Turbulent Flows”, AIAA Journal, 37, 2, (1999) 185-191. 

[18] Burdyshaw, C. E., Anderson, W. K., “A General and Extensible Unstructured Mesh 
Adjoint Method”, AIAA J. of Aerospace, Computing, Information, and Communication, 2, 
10,401-413 (2005). 

[19] Burdyshaw, C. E., “Achieving Automatic Concurrency Between Computational Field 
Solvers and Adjoint Sensitivity Codes,” Ph.D. Thesis, University of Tennessee, 
Chattanooga, May, 2006. 

[20] Jameson, A., “Aerodynamic Design Via Control Theory”, J. Scientific Computing, 3, 
(1998)233-260. 

[21] Jameson, A., Alonso, J. J., Reuther, J., Martinelli, L., Vassberg, J. C., “Aerodynamic 
Shape Optimization Techniques Based on Control Theory”, AIAA Paper No. 98-2538, 
(1998). 

[22] Mohammadi, B., “Optimal Shape Design, Reverse Mode of Automatic Differentiation 
and Turbulence”, AIAA Paper No. 97-0099, (1997). 


37 


[23] Nielsen, E. J., Anderson, W. K., “Recent Improvements in Aerodynamic Optimization 
on Unstructured Meshes”, AIAA J., 40, 6, (2002) 1155-1163. 

[24] Nielsen, E. J., Anderson, W. K., “Aerodynamic Design Optimization On Unstructured 
Meshes Using the Navier-Stokes Equations”, AIAA J., 37, 11, (1999) 1411-1419. 

[25] Nielsen, E. J., “Aerodynamic Design Sensitivities on an Unstructured Mesh Using the 
Navier-Stokes Equations and a Discrete Adjoint Formulation”, Ph.D. Dissertation, Virginia 
Polytechnic Institute and State University, (1998). 

[26] Nielsen, E. J., Kleb, W. L., “Efficient Construction of Discrete Adjoint Operators on 
Unstructured Grids by Using Complex Variables”, AIAA J., 44, 4 (2005) 827-836. 

[27] Park, M., “Adjoint-Based, Three-Dimensional Error Prediction and Grid Adaptation”, 
AIAA Paper No. 2002-3286, (2002). 

[28] Kapadia S., Anderson W. K., Elliott L., Burdyshaw C., “Adjoint method for solid- 
oxide fuel cell simulations”, J. Power Sources, 166 (2007) 376-385. 

[29] Kapadia S., Anderson W. K., Elliott L., Burdyshaw C., “Adjoint based Sensitivity 
Analysis and Error Correction Methods applied to Solid Oxide Fuel Cells”, ASME J. Fuel 
Cell Sci. Tech., 6, 2, (2009). 

[30] Kapadia S., Anderson W. K., “Sensitivity Analysis for Solid Oxide Fuel cells using a 
Three-Dimensional Numerical Mode”, J. Power Sources, 189 (2009) 1074-1082. 

[31] Kapadia S., Anderson W. K., Burdyshaw, C., “Channel shape optimization of solid 
oxide fuel cells using advanced numerical techniques”, Computers and Fluids, 41,1 (2011) 
41-50. 


[32] Lin C., Huang L., Chiang L., Chyou Y., “Effects of clamping load on the thermal 
stress distribution in a planar SOFC with compressive sealing”, ECS Transactions, 25, 2 
(2009) 349-358. 

[33] Lin C., Huang L, Chiang L., Chyou Y., “Thermal stress analysis of planar solid oxide 
fuel cell stacks: Effects of sealing design”, J. Power Sources 192 (2009) 515-524. 

[34] Chiang L, Liu H., Shiu Y., Lee C., Lee R., “Thermo-electrochemical and thermal 
stress analysis for an anode-supported SOFC cell”, Renewable Energy 33 (2008) 2580- 
2588. 

[35] Iqbal G., Pakalapati S., Blancas F., Guo H., Celik I., Kang B., “PEN structure thermal 
stress analysis for planar-SOFC configurations under practical temperature field”. 


38 


Advances in Materials Science for Environmental and Nuclear Technology II, Ceramic 
Transactions 227 (2011) 61-68. 

[36] Jiang T., Chen T., “Thermal-stress analyses of an operating planar solid oxide fuel 
cell with the bonded compliant seal design”, Int. J. Hydrogen Energy 34 (2009) 8223- 
8234. 

[37] Weil K., Koeppel B., “Comparative finite element analysis of the stress-strain states in 
three different bonded solid oxide fuel cell seal designs”, J. Power Sources 180 (2008) 
343-353. 

[38] Chiang L., Liu H, Shiu Y., Lee C., Lee R., “Thermal stress and thermo¬ 
electrochemical analysis of planar anode-supported solid oxide fuel cell: Effect of anode 
porosity”, J. Power Sources 195 (2010) 1895-1904. 

[39] ABAQUS webpage: http://www.3ds.com/products/simulia/portfolio/abaqus/abaqus- 
portfolio/ 

[40] FLUENT webpage: 

http://www.ansys.com/Products/Simulation+Technology/Fluid+Dynamics/Fluid+Dynamic 
s+Products/AN S YS+Fluent 

[41] ANSYS webpage: http://www.ansys.com 

[42] STAR-CD webpage: http://www.cd-adapco.com/products/star_cd/ 

[43] MARC webpage: http://www.mscsoftware.com/product/marc 

[44] Wang, Y., Yoshiba, F., Watanbe, T., Weng, S., “Numerical analysis of 
electrochemical characteristics and heat/species transport for planar porous-electrode- 
supported SOFC”, J. Power Sources, 170, (2007) 101-110. 

[45] Hwang, J. J., Chen, C. K., Lai, D. Y., “Computational analysis of species transport 
and electrochemical characteristics of a MOLB-type SOFC”, J. Power. Sources, 140, 
(2005) 235-242. 

[46] Noren, D. A., Hoffman, M. A., “Clarifying the Butler-Volmer equation and related 
approximation for calculating activation losses in solid oxide fuel cell models”, J. Power 
Sources, 152, (2005) 175-181. 

[47] Costamagna, P., Selimovic, A., Borghi, M., Agnew, G., “Electrochemical model of 
the integrated planar solid oxide fuel cell (IP-SOFC)”, Chem. Eng. J., 102, (2004) 61-69. 


39 


[48] Saad, Y., Schultz, M. H., “GMRES: A Generalized Minimal Residual Algorithm for 
Solving Nonsymmetric Linear Systems”, SIAM J. Sci. Stat. Comput., 7, (1986) 856-869 

[49] Grama A., Gupta, A., Karypis, G., Kumar, V., “Introduction to Parallel Computing”, 
Second Edition, Addison Wesley, Jan. 16, 2003. 

[50] METIS webpage: http://glaros.dtc.umn.edu/gkhome/views/metis 

[51] POINTWISE webpage: http://www.pointwise.com 

[52] Roe, P. L., “Characteristic-based schemes for the Euler equations”, Ann. Rev. Fluid 
Mech., 18,(1986) 337-365. 

[53] Busby M., “Steps Toward More Accurate and Efficient Simulations of Reactive 
Flows”, Ph.D. Thesis, Mississippi State University, August 1997. 


40 


